rm(list=ls())
gc()

library(dplyr)
library(RSwissMaps)
library(ggpubr)

####### Plot the results #####


####### Party competition space ######

###### Simple model 
load("model/Model party ideal position.rda")


thetas <- extract(model1, pars = "theta")
thetas <- as.data.frame(thetas)

png("Plots/Figure 1.png", width=2000, height=1500, res=300)
ggplot(thetas[1:5000,])+ ### Plot only draws from the first chain
  geom_density(aes(x=theta.1, color = "CVP", fill ="CVP"), alpha = .3)+
  geom_density(aes(x=theta.2, color = "FDP", fill = "FDP"), alpha = .3)+
  geom_density(aes(x=theta.3, color = "GPS", fill = "GPS"), alpha = .3)+
  geom_density(aes(x=theta.4, color = "SPS", fill = "SPS"), alpha = .3)+
  geom_density(aes(x=theta.5, color = "SVP", fill = "SVP"), alpha = .3)+
  scale_color_manual(name = "", values = c("orange", "blue", "green", "red", "dark green"))+
  scale_fill_manual(name = "", values = c("orange", "blue", "green", "red", "dark green"))+
  theme_minimal()+
  xlab(expression(theta[j]))+
  ylab("Density")+
  theme(legend.position = "bottom")
dev.off()


####### Canton and municipalities in party competition space
load("data/Results position municipalities and cantons.Rda")


#### Figure 2
data_canton_results <- data_results_canton_mun$data_canton_results
data_mun_results <- data_results_canton_mun$data_mun_results

data_canton_results$name <- data_canton_results$number_mun
data_mun_results$bfs_nr <- as.numeric(data_mun_results$number_mun)

### Dynamic models 
data_canton_results_time <- data_results_canton_mun$data_canton_results_time
data_mun_results_time <- data_results_canton_mun$data_mun_results_time

data_canton_results_time$name <- data_canton_results_time$number_mun
data_mun_results_time$bfs_nr <- as.numeric(data_mun_results_time$number_mun)


##### Load the template
can_plot <- can.template(2016)
mun_plot <- mun.template(2016)

can_plot1 <-  merge(can_plot, data_canton_results, by = "name")

mun_plot1 <-  merge(mun_plot, data_mun_results, by = "bfs_nr")

png("Plots/Figure 2.png", width=3000, height=1500, res=300)
ggarrange(
can.plot(bfs_id = can_plot1$bfs_nr, data = can_plot1$rand, year = 2016, boundaries_size = 0.1,
         legend_title = expression(u["1g"]), boundaries_color = "white")+
  scale_fill_gradient2(
    name = expression(u["1g"]),
    low = "red",
    mid = "light grey",
    high = "dark green",
    midpoint = 0,
    space = "Lab",
    na.value = "white",
    guide = "colourbar",
    aesthetics = "fill"
  ),

mun.plot(bfs_id = mun_plot1$bfs_nr, data = mun_plot1$rand, year = 2016, boundaries_size = 0.1,
         legend_title = expression(u["1g"]), boundaries_color = "white")+
  scale_fill_gradient2(
    name = expression(u["1g"]),
    low = "red",
    mid = "light grey",
    high = "dark green",
    midpoint = 0,
    space = "Lab",
    na.value = "white",
    guide = "colourbar",
    aesthetics = "fill"
  ), ncol=2)
dev.off()


#### Figure 3

order_can_plot <- order(can_plot1$rand)

can_plot1$order <- NA

for (i in 1:nrow(can_plot1)) {
  can_plot1[i,]$order <- which(order_can_plot==i)
}

can_plot1 <- can_plot1 %>% 
  mutate(col = ifelse(rand<0 & high<0, "Significantly more progressive \n than the Swiss population",
                      ifelse(rand<0 & high>0, "More progressive than the Swiss population \n not significant", 
                             ifelse(rand>0 & low>0, "Significantly more conservative \n than the Swiss population", 
                                    ifelse(rand>0 & low<0, "More conservative than the Swiss population \n not significant", NA)))))

order_mun_plot <- order(mun_plot1$rand)

mun_plot1$order <- NA

for (i in 1:nrow(mun_plot1)) {
  mun_plot1[i,]$order <- which(order_mun_plot==i)
}

mun_plot1 <- mun_plot1 %>% 
  mutate(col = ifelse(rand<0 & high<0, "Significantly more progressive \n than the Swiss population",
                      ifelse(rand<0 & high>0, "More progressive than the Swiss population \n not significant", 
                             ifelse(rand>0 & low>0, "Significantly more conservative \n than the Swiss population", 
                                    ifelse(rand>0 & low<0, "More conservative than the Swiss population \n not significant", NA)))))


png("Plots/Figure 3.png", width=3000, height=3000, res = 300)
ggarrange(
ggplot(can_plot1)+
  geom_point(aes(x=rand, y = order, color = col))+
  geom_errorbarh(aes(xmin = low, xmax=high, y = order, color = col), alpha = .2)+
  scale_color_manual(values = c("lightgreen", "red", "dark green", "dark red"))+
  theme_minimal()+
  xlab(expression(u["1g"]))+
  ylab("Cantons order (progressive to conservative)")+
  theme(legend.position = "bottom", legend.title = element_blank()), 


can.plot(bfs_id = can_plot1$bfs_nr, data= as.factor(can_plot1$col), 2016, continuous = FALSE, color_discrete = c("qual", "2"),
         boundaries_size = 0.1)+
  scale_fill_manual(values = c("lightgreen", "red", "darkgreen", "darkred"), na.value = "white",na.translate = FALSE)+
  theme(legend.title = element_blank()),



ggplot(mun_plot1)+
  geom_point(aes(x=rand, y = order, color = col))+
  geom_errorbarh(aes(xmin = low, xmax=high, y = order, color = col), alpha = .2)+
  scale_color_manual(values = c("lightgreen", "red", "darkgreen", "darkred"))+
  theme_minimal()+
  xlab(expression(u["1g"]))+
  ylab("Municipalities order (progressive to conservative)")+
  theme(legend.position = "bottom", legend.title = element_blank()),

mun.plot(bfs_id = mun_plot1$bfs_nr, data= as.factor(mun_plot1$col), 2016, continuous = FALSE, color_discrete = c("qual", "2"),
         boundaries_size = 0.1)+
  scale_fill_manual(values = c("lightgreen", "red", "darkgreen", "darkred"), na.value = "white",na.translate = FALSE)+
  theme(legend.title = element_blank()), nrow=2, ncol=2, common.legend = T, legend = "bottom")
dev.off()


######## Figure 4 #######

###### Dynamic model
load("model/Model party ideal position_by legislature.rda")

thetas <- summary(model2, pars = "theta2")
thetas <- as.data.frame(thetas)

thetas$j <- rep(c("CVP", "FDP", "GPS", "SPS", "SVP"), 16)
thetas$t <- rep(1:16, each=5)



png("Plots/Figure 4.png", width=2000, height=3000, res=300)
ggarrange(
  ggplot(thetas)+
    geom_line(aes(x=t+35, y=c_summary.50..chain.1, color = j))+
    geom_line(aes(x=t+35, y=c_summary.50..chain.1 - 1.96 * c_summary.sd.chain.1, color = j), linewidth=.2)+
    geom_line(aes(x=t+35, y=c_summary.50..chain.1 + 1.96 * c_summary.sd.chain.1, color = j), linewidth=.2)+
    theme_minimal()+
    scale_color_manual(values = c("orange", "blue", "green", "red", "darkgreen"))+
    xlab("Legislature")+
    ylab(expression(theta[jt]))+
    theme(legend.position = "bottom", legend.title = element_blank()),
ggplot(data_canton_results_time)+
  geom_line(aes(x=as.numeric(t), y=rand, color = name), alpha = 0.5)+
  scale_color_manual(values = rep("grey", length(unique(data_canton_results_time$name))))+
  theme_minimal()+
  xlab("Legislature")+
  labs(subtitle = "Cantons")+
  ylab(expression(u[gt(i)]))+
  theme(legend.position = "none"),

ggplot(data_mun_results_time)+
  geom_line(aes(x=as.numeric(t), y=rand, color = number_mun), alpha=0.3)+
  scale_color_manual(values = rep("grey", length(unique(data_mun_results_time$number_mun))))+
  theme_minimal()+
  labs(subtitle = "Municipalities")+
  xlab("Legislature")+
  ylab(expression(u[gt(i)]))+
  theme(legend.position = "none"), align = "hv", nrow=3)
dev.off()


### Figure 5 ###


rm(list=ls())
gc()
load("model/Model IRT with parties and cantons.Rda")
load("model/Model IRT with parties and municipalities.Rda")

betas_party_canton <- summary(model_party_canton, pars = "beta")
betas_party_canton <- as.data.frame(betas_party_canton)

thetas_party_canton <- summary(model_party_canton, pars = "theta")
thetas_party_canton <- as.data.frame(thetas_party_canton)

betas_party_mun <- summary(model_party_mun, pars = "beta")
betas_party_mun <- as.data.frame(betas_party_mun)

thetas_party_mun <- summary(model_party_mun, pars = "theta")
thetas_party_mun <- as.data.frame(thetas_party_mun)


thetas_party_canton$j <- 1:31
thetas_party_mun$j <- 1:2177

thetas_party1 <- thetas_party_canton[thetas_party_canton$j>26,]
thetas_party1$j <- thetas_party1$j - 26

thetas_party2 <- thetas_party_mun[thetas_party_mun$j>2172,]
thetas_party2$j <- thetas_party1$j - 2172

thetas_canton <- thetas_party_canton[thetas_party_canton$j<27,]
thetas_mun <- thetas_party_mun[thetas_party_mun$j<2173,]


######### Extract parameters from the model with party #####
load("model/Model party ideal position.rda")
# load("model/Model party ideal position_by legislature.rda")


betas_o <- summary(model1, pars="beta")
betas_o <- as.data.frame(betas_o)

thetas_o <- summary(model1, pars="theta")
thetas_o <- as.data.frame(thetas_o)



#### Compare discrimination parameter from the two models 
data_compare <- as.data.frame(list(original = betas_o$c_summary.50..chain.1,
                                   party_mun = betas_party_mun$c_summary.50..chain.1))

data_compare <- data_compare %>% 
  mutate(diff = abs(original - party_mun))

##### Static 
png("Plots/Figure 5.png", width=3000, height=1500, res=300)
ggarrange(
  ggplot()+
    geom_point(aes(x=thetas_o$c_summary.50..chain.1, y=thetas_party1$c_summary.50..chain.1))+
    geom_errorbar(aes(x=thetas_o$c_summary.50..chain.1, 
                      ymin = thetas_party1$c_summary.50..chain.1 - 1.96 * thetas_party1$c_summary.sd.chain.1,
                      ymax = thetas_party1$c_summary.50..chain.1 + 1.96 * thetas_party1$c_summary.sd.chain.1), alpha = .2)+
    geom_errorbarh(aes(y=thetas_party1$c_summary.50..chain.1, 
                       xmin = thetas_o$c_summary.50..chain.1 - 1.96 * thetas_o$c_summary.sd.chain.1,
                       xmax = thetas_o$c_summary.50..chain.1 + 1.96 * thetas_o$c_summary.sd.chain.1), alpha = .2)+
    xlab("Parties' positions in the party competition space")+
    ylab("Parties' positions in the party canton space")+
    theme_minimal(),
  
  
  
  ggplot()+
    geom_point(aes(x=thetas_o$c_summary.50..chain.1, y=thetas_party2$c_summary.50..chain.1))+
    geom_errorbar(aes(x=thetas_o$c_summary.50..chain.1, 
                      ymin = thetas_party2$c_summary.50..chain.1 - 1.96 * thetas_party2$c_summary.sd.chain.1,
                      ymax = thetas_party2$c_summary.50..chain.1 + 1.96 * thetas_party2$c_summary.sd.chain.1), alpha = .2)+
    geom_errorbarh(aes(y=thetas_party2$c_summary.50..chain.1, 
                       xmin = thetas_o$c_summary.50..chain.1 - 1.96 * thetas_o$c_summary.sd.chain.1,
                       xmax = thetas_o$c_summary.50..chain.1 + 1.96 * thetas_o$c_summary.sd.chain.1), alpha = .2)+
    xlab("Parties' positions in the party competition space")+
    ylab("Parties' positions in the party municipality space")+
    theme_minimal(), ncol = 2)
dev.off()


######### Figure 6 Compare item parameter value ######

png("Plots/Figure 6.png", width=3000, height=1500, res=300)
ggarrange(
  ggplot()+
    geom_point(aes(y=betas_party_canton$c_summary.50..chain.1, x=betas_o$c_summary.50..chain.1))+
    geom_errorbarh(aes(xmin= betas_o$c_summary.50..chain.1-1.96 * betas_o$c_summary.sd.chain.1,
                       xmax=betas_o$c_summary.50..chain.1+1.96 * betas_o$c_summary.sd.chain.1,
                       y=betas_party_canton$c_summary.50..chain.1), alpha = .2)+
    geom_errorbar(aes(ymin= betas_party_canton$c_summary.50..chain.1-1.96 * betas_party_canton$c_summary.sd.chain.1,
                      ymax=betas_party_canton$c_summary.50..chain.1+1.96 * betas_party_canton$c_summary.sd.chain.1,
                      x=betas_o$c_summary.50..chain.1), alpha = .2)+
    xlab("Discrimination of ballots in the party competition space")+
    ylab("Discrimination of ballots in the party canton space")+
    theme_minimal(),
  
  ggplot()+
    geom_point(aes(y=betas_party_mun$c_summary.50..chain.1, x=betas_o$c_summary.50..chain.1))+
    geom_errorbarh(aes(xmin= betas_o$c_summary.50..chain.1-1.96 * betas_o$c_summary.sd.chain.1,
                       xmax=betas_o$c_summary.50..chain.1+1.96 * betas_o$c_summary.sd.chain.1,
                       y=betas_party_mun$c_summary.50..chain.1), alpha = .2)+
    geom_errorbar(aes(ymin= betas_party_mun$c_summary.50..chain.1-1.96 * betas_party_mun$c_summary.sd.chain.1,
                      ymax=betas_party_mun$c_summary.50..chain.1+1.96 * betas_party_mun$c_summary.sd.chain.1,
                      x=betas_o$c_summary.50..chain.1), alpha = .2)+
    xlab("Discrimination of ballots in the party competition space")+
    ylab("Discrimination of ballots in the party municipality space")+
    theme_minimal(), ncol = 2, align = "hv")
dev.off()  

